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ABSTRACT 

We study the interaction between massive planets and a gas disc with a mass in the 
range expected for protoplanetary discs. We use SPH simulations to study the orbital 
evolution of a massive planet as well as the dynamical response of the disc for planet 
masses between 1 and 6 Mj and the full range of initial relative orbital inclinations. 

We find that gap formation can occur for planets in inclined orbits as well as 
for coplanar orbits as expected. For given planet mass, a threshold relative orbital 
inclination exists under which a gap forms. This threshold increases with planet mass. 
Orbital migration manifest through a decreasing semi-major axis is seen in all cases. 

At high relative inclinations, the inclination decay rate increases for increasing 
planet mass and decreasing initial relative inclination as is expected from estimates 
based on dynamical friction between planet and disc. For an initial semi-major axis 
of 5 AU and relative inclination of io — 80°, the times required for the inclination to 
decay by 10° is ^ 10^ yr and ^ 10^ yr for 1 Mj and 6 Mj respectively, these times 
scaling in the usual way for larger initial orbits. For retrograde planets, the inclination 
always evolves towards coplanarity with the disc, with the rate of evolution being 
fastest for orbits with iq — > 180°. The indication is thus that, without taking account 
of subsequent operation of phenomena such as the Lidov-Kozai effect, planets with 
mass 1 Mj initiated in circular orbits with semi-major axis ^ 5 AU and iq ~ 90° might 
only just become coplanar, as a result of frictional effects, within the disc lifetime. In 
other cases highly inclined orbits will survive only if they are formed after the disc 
has mostly dispersed. 

Planets on inclined orbits warp the disc by an extent that is negligible for 1 Mj 
but increases with increasing mass becoming quite significant for a planet of mass 
6 Mj. In that case, the disc can gain a total inclination of up to 15° together with a 
warped inner structure with an inclination of up to ^ 20° relative to the outer part. 
We also find a solid body precession of both the total disc angular momentum vector 
and the planet orbital momentum vector about the total angular momentum vector, 
with the angular velocity of precession decreasing with increasing relative inclination 
as expected in that case. 

Our results illustrate that the influence of an inclined massive planet on a pro- 
toplanetary disc can lead to significant changes of the disc structure and orientation 
which can in turn affect the orbital evolution of the planet significantly. A three- 
dimensional treatment of the disc is then essential in order to capture all relevant 
dynamical processes in the composite system. 
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1 INTRODUCTION 

In the past few decades, the number of detected extraso- 
lar planets around other main sequence stars has increased 
dramatically, so that at the time of writing, more than 840 
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planets have been discovered. W ell studied planet formation 
scenarios such as core accretion (lMizunolll98(]|:|Pollack et all 
Il996h or disc fragmentation (jMaver et al.ll2002l " involve the 
planet forming in a disc with the natural expectation that 
the orbit should be coplanar. However, Rossiter McLaugh- 
lin measurements of close orbiting hot Jupiters indicate that 
around 40% have angular momentum vectors significantly 
misaligne d with the angular velocity vector of th e central 
star (e.g. iTriaud et al.ll201(]| : lAlbrecht et al.ll2012l ). As the 
stellar and disc angular velocities are naturally expected to 
be aligned, this would imply an inclination of the planet or- 
bit r elative to the na scent protoplanetary disc. But note that 
fe.g. lLai et aLll201lh have proposed that a magnetic interac- 
tion between the central star and protoplanetary disc could 
lead to the stellar angular momentum vector being mis- 
aligned with that of the disc, in which case misali gned plan- 
etary orbits need not be inferred. However, Watson et al.l 
l|201ll ) have compared stellar rotation axis inclination angles 
with measured debri-disc inclinations and found no evidence 
for misalignment. 

Several scenarios have been proposed to explain the 
origin of misalignments between the planetary orbital an- 
gular momentum vector and the stellar rotation axis for 
close in planets. The first involves excitation of very 
high eccentricities, either through the Lidov-Kozai ef- 
fect indu ced by the in teraction with a distant com pan- 
ion (e.g. iFabrvckv fc"Tremaine ,2007]; IWu et al.l [20071 ) . or 
thro ugh planet-planet scattering or chaotic interactions 
fe.g. IWeidenschillin g fc Marzari' 199 61: iRasio fc Fordlll996l : 
iPapaloizou fc Terq ucm 2001; Nagasa wa et al.ll2008l ~This is 
is then followed by orbital circularization due to interaction 
with the central star, alth ough aspects of this process are not 
yet understood in detail (jPawson et al.ll2012l ). It is possible 
that high eccentricities and inclinations for massive planets 
are generated when th e disc is presen t . It is then important, 
as has been noted by IPawson et al.l ()2012[ ). to understand 
the interactions with the disc and consequent time scales for 
realignment. 

Another meth od for pr oducing high inclinations has 
been indicated bv iThomme s & Lis sauer. ( 200 3), who have 
studied the evolution of two giant planets in resonance, with 
an approximate analytical expression for the influence of a 
gas disc in producing orbital migration. Their calculations 
suggest that resonant inclination excitation can occur when 
through the interaction between the planets the eccentricity 
of the inner planet reaches a threshold value of ei 0.6. In- 
clinations gained by the resonant pair of planets can reach 
values up to ~ 60°. The effects of disc planet interactions 
may be significant here also. 

Another explanation for the origin of misaligned plan- 
ets is connected with the possibility that the orientation 
of the disc changes, possibly due to lack of alignment of 
the source material or due to gravitational encounters with 
passing stars, ei t her before or afte r the planet forms (e.g. 
iBate et al.l I2OI0I : iThies et al1l201ll '). However, up to now, 
there is no observational evidence for such disc misalign- 
ment (|Watson et al.ll201ll ). 

While the formation of planets on misaligned orbits 
has been a frequent topic of debate, the evolution of mas- 
sive planets on misaligned orbits interacting with a disc is 
still not well understood. Many hydrodynamical simulations 
have been performed in order to study the influence of a pro- 



toplanetary disc on the evolution of planets. The majority of 
these simulations have considered planets that are coplanar 
with the disc and been two di mensional and used grid-based 
methods (see IPapaloizou et aL 2007 ; ,Klev fc Nel son 20l"3. 
for a review and references therein). 

Planets with initial orbital inc linations with respect 
to a disc have bee n studied (e.g. ICresswelT et all l2007l : 
iBitsch fc Kievll201ll ). In these works the evolution of plan- 
ets with masses up to a maximum of IMj with different 
initial eccentricities and relatively small initial inclinations 
up to a maximum of 15° interacting with three-dimensional 
isothermal and radiative discs are considered. They showed 
damping for both eccentricity and inclination with the plan- 
ets circularising in the disc after a few hundred orbits. 

In this paper we are interested in extending studies 
of disc planet interactions involving planets with orbital 
planes misaligned with that of the disc to consider larger 
planet masses of up to several Jupiter masses and the full 
range of inclinations. In such cases the interaction with the 
disc is more likely to be described by dynamical friction 
than being the result of the application of resonant torques 
that is applicable in the coplanar case (see for example 
Papaloizou fc Larwoodll200(]| ; lReiDll2012l ; iTevssandier et all 



20131) . Numerical simulations of wa rped discs in c l ose b i- 
nary systems were perfor med by e.g. lLarwood et al] l|l996l ); 
iFragncr fc NelsonI ([201(1). These works indicate that thick 
discs with low viscosities have efflcient warp communica- 
tion which allows them to precess as rigid bodies with little 
warping or twisting while thinner discs can develop twists 
before reaching a state of rigid-body precession. Possible 
warping of a disc by giant planets has yet to be studied 
in detail. It has been est imated from linearized calculations 
(|Tevssandier et ai]|2013l ) that disc warping should be a small 
effect for sub Jovian mass planets, this being confirmed by 
the results of this paper. However, it turns out to be signif- 
icant for larger planet masses that can produce significant 
nonlinear effects, therefore we study these effects in this pa- 
per. 

It has been shown that Smoothed particle hydrody- 
namics (SPH) simulations are capable of bei ng applied to 
the problem of planet-disc int eraction (e.g. I S chafer et al.l 
|2004 ; Ide Val-Borro et all bOOel ) though in comparison to 
grid-based simulation methods, they incur greater compu- 
tational expense. But in contrast to grid-based methods, as 
they adopt a Lagrangian approach, SPH simulations can 
be readily used to simulate a gas disc with a free bound- 
ary in three dimensions with the disc having the freedom to 
change its shape at will, thu s we a dopt this approach here. 
Similarly, iMarzari fc NelsonI (|2009l ) studied the interaction 
between a Jupiter mass planet and a disc using the SPH 
method. Initial eccentricities ranging from to 0.4 and an 
initial inclination of 20° were considered. The disc is simu- 
lated by the SPH code VINE that is similar to GADGET-2 
and a locally isothermal equation of state. They find that 
both the eccentricity and the inclination are damped on a 
timescale much lower that the disc lifetime (~ 10® yr). 

In this paper we consider a range of planet masses and 
the full range of orbital inclinations. In addition to studying 
the warping response of the disc, we estimate the timescales 
for orbital evolution and attainment of coplanarity. We also 
make a preliminary study of the potential exchange of in- 
clination and eccentricity that might be expected for high 
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inclination orbits through the operation of an adaption of 
the Lidov-Kozai mechanism. The plan of the paper is as 
follows: 

In Section [5J we describe our simulation code with a 
description of the modified locally isothermal equation of 
state we adopted fSection l2.1[l . Some details concerning the 
smoothing length and artificial viscosity are given in Sec- 
tion [2]2] In Section [3l we describe the general setup for the 
disc, planet and the central star. In Section 14] we give a 
brief overview of the interaction between a massive planet 
and disc. We give an analytic estimate of the evolution time 
scale for a planet interacting with a disc based on consider- 
ations of dynamical friction, in order to compare with our 
simulations, in Section[4]2] In addition, we consider the pos- 
sibility of the Lidov-Kozai effect (Section 14. 3p . In Section O 
we discuss numerical results for a massive planet on a copla- 
nar circular orbit. We study the evolution of planets initi- 
ated on inclined orbits in Section [5] considering the evolution 
of both initially prograde and retrograde circular orbits. In 
Section [T] we study disc warping and precession. In Section 
[S] we study a planet starting on an eccentric inclined orbit 
finding evidence for some exchange of inclination and eccen- 
tricity as might be expected from a Lidov-Kozai like process. 
The behaviour of an inclined planet interacting with discs 
of different masses is considered in Section fS.ll Finally the 
results are summarised and discussed in Section (9) 



2 SIMULATION DETAILS 

We have performed simulations using a mo dified version o f 
the publically available code GADGET-2 (|Springelll2005[ ). 
GADGET-2 is a hybrid N-body/SPH code capable of mod- 
elling both fluid and distinct fixed or orbiting massive bod- 
ies. In our case the central star, of mass M* , is fixed and the 
planet, of mass Mp, orbits as a distinct massive body. We 
adopt spherical polar coordinates (r, 6, (j>) with origin at the 
centre of mass of the central star. 

The gaseous disc is represented by SPH particles. An impor- 
tant issue in N-body/SPH simulations is the choice of the 
gravitational softening lengths. Ideally, gravitational force 
computations without softening should be used to deter- 
mine the gravitational interactions between particles. How- 
ever, for numerical reasons, a softening length has to be 
introduced. This is because SPH particles would experience 
very high gravitational accelerations in the vicinity of mas- 
sive objects. In turn, this would result in unacceptably small 
time steps that prevent simulations from effectively proceed- 
ing. For this reason, we adopt gravitational softening lengths 
that are sufficiently large to enable the calculation to pro- 
ceed, and in the case of the planet, smaller than physically 
relevant length scales such as the disc height H and the plan- 
etary Hill radius. It can be associated with a physical size 
of material considered to be bound to the planet. Here we 
study planets of masses Mp ^ 1 Mj with an orbital length 
scale of Rp = 5 AU and discs with aspect ratio H/r — 0.05. 
Thus the HiU radius defined as Rh = Rp (Mp/(3M.))^/^ is 
0.35 AU for a Jupitermass planet at 5 AU. At r = 5 AU, 
H — 0.25 AU. In our simulations, the central star and 
the planet are unsoftened when interacting with each other 
allowing them to undergo close encounters accurately, al- 
though they do not actually occur in our simulations. For 



the computation of the gravitational interaction between the 
massive bodies and the SPH particles, we adopted fixed soft- 
ening lengths e» = 0.4 AU and Sp — 0.1 AU for the central 
star and the planet, respectively. The former value controls 
conditions near the central star. These are not significant 
for the dynamics of interest in our simulations which only 
depend on the value of the central stellar mass. The latter 
value Ep controls conditions close to the planet and is signif- 
icantly smaller than both Rh{1 Mj, 5 AU) and H{5 AU). 

The star and the planet are allowed to accrete gas par- 
ticles that approach them very closely. The reason for doing 
this is to prevent the possibility of strong unresolved forces 
occuring close to massive objects. In this respect it works 
in tandem with softening. In order to impleme nt accretion, 
we followed the procedure of lBate et al. I (fl995l ') for so-called 
"sink particles". This was applied such that for the star and 
planet, the outer accretion radii were fixed during the simu- 
lation to be Racer,- ~ 7 X 10^^ cm and Raccr,p = 7 x lO^'' cm 
respectively. SPH particles that come within the accretion 
radius of a sink particle, are accreted if they fuUfiU all con- 
ditions for accretion. Firstly, the particle must be bound to 
the sink particle. Secondly, its specific angular momentum 
about the sink particle must be less than that required for it 
to form a circular orbit at Race about the sink particle. These 
conditions ensure that particles which would normally leave 
the accretion radius are not accreted. Thirdly, the particle 
must be more tightly bound to the candidate sink particle 
than to any other sink particle. In our case this is to allow 
for the possibility that, given the very different masses of 
the star and planet, a SPH particle could pass through the 
accretion radius of the planet and yet be accreted by the 
central star. 

An inner accretion radius can also be define d such that 
SPH particles are accreted regardless of the tests (|Bate et al.l 
1995). By using an inner accretion radius, particles are pre- 
vented from being accelerated by very close encounters with 
the sink particle tha t might otherwise dramatically slow 
down the simulation. iBate et"ai] (|l995l ) suggests that the 
inner accretion radius should be 10-100 times smaller than 
the outer accretion radius. Due to our very small outer ac- 
cretion radii, the choice of the inner accretion radii to be 0.5 
Race was found to be adequate for this purpose. 

In practice, for the equations of state and softening pa- 
rameters we used for most of our simulations, because of 
the relatively small accretion radii adopted, the accretion of 
gas particles was found to play only a minor role, producing 
negligible changes to the masses of the planet and star (see 
also the discussion in section [STTJ . 



2.1 Equation of state 

For the hydrodynamical simulations, we applied a locally 
iso thermal equation of sta te with the modification described 
by iPeplinski et al' (2008*). The reason for modification is 
because, as Pcplinski ct al. ( 200& ) have pointed out, if the 
temperature characteristic of the main disc is adopted in 
the vicinity of the planets, a high concentration of gas lead- 
ing to rapid accretion develops in the close vicinity of mas- 
sive planets. Accordingly, we adopt a sound speed given by 
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l|Peplinski et al.ll2008l ) 

hsVs hpVp 



[{hsTs 



Here rs 



— Tsl and r 



(1) 



Vnl are the distances to 



the central star and the planet respectively, and ils and Q,p 
are the angular velocities in the circumstellar and circum- 
planetary discs. For pure Keplerian motion they are given 
by 



GAL 



and ilr, 



GMp 



respectively. 



(2) 



The disc aspect ratio is hs — H/vs with H being the disc 
scale height. The parameter n was taken to be n = 3.5. 
For the simulations we take H/rs — Cs/v^ — 0.05, where 
Vip = rsQs is the rotational velocity. As rp — >■ 0, the sound 
speed in the circumplanetary disc is hpVpQp and accord- 
ingly hp is the aspect ratio in that limit. This should be 
chosen large enough so that pressure forces prevent exces- 
sive gas accumulation on the planet. At the same time, it 
must be chosen large enough to prevent values of the sound 
speed that are smaller than the unmodified value, at loca- 
tions that are a mod erate distance away from the planet. 
IPeplinski et al.l (|2008l ^ suggest hp ^ 0.4. A discussion of the 
effects resulting from different choices of hp can be found in 
appendix|Al For our simulations, we adopted hp = 0.6 which 
was large enough to prevent unrealistically fast gas accretion 
while not causing the sound speed to become smaller than 
the unmodified value in all cases. 



2.2 Smoothing length and artificial viscosity 

For our SPH calculations, the smoothing length was ad- 
justed so that the number of nearest neighbours to any 
particle contained within a sphere of radius equal to the 
local smoothing length was 40 ± 5. The pressure is given 
by p = pCg. Thus apart from the vicinity of the planet, the 
temperature in the disc is oc r^^. The artificial viscosity 
parameter a o f GADGET-2 (see equations (9) and (14) of 
ISpringelllioOsI ) was taken to be a = 0.5. 

The artificial viscosity is modified by the application 
of a viscosity-limiter to reduce artificially induced angular 
momentum transport in the presence of shear fiows. This 
is especially important for the study of Keplerian discs. In 
order to calibrate the diffusive effects resulting from this 
viscosity, we applied the ring spreading test to determine 
an effective kinematic viscosity meas ured through the ass 
parameter (jShakura fc SunvaevlllQT^ ). Details are given in 
appendix [B] A value of a = 0.5, showed the best match to 
the analytic ring spreading solution with ass ~ 0.02 and 
was hence used in all simulations. SPIf simulations accord- 
ingly model a viscous disc, with a viscosity that behaves 
like a conventional Navier Stokes viscosity in this context. 
It should be noted that the effective viscosity in protoplan- 
etary discs is likely to involve magnetic fields and behave 
differently in some contexts. However many of the simula- 
tions undertaken here involve impulsive high velocity inter- 
actions between planet and disc that are characteristically 
described by dynamical friction (see below) and for which 
the disc viscosity is not expected to play a major role. The 
remaining simulations involve coplanar, inwardly migrating. 



gap forming planets for which the interaction with the disc 
is nonlinear and long range. The characteristic behaviour 
again is not e xpected to be sensitive to det ails of the viscos- 
ity (e. g. iNelson &: Papaloizou. 2003l : .Baruteau et al.ll2011i V 



3 INITIAL CONDITIONS 

We study a system composed of a central star of one solar 
mass M0, a gaseous disc and a massive planet. The disc 
is set up such that the angular momentum vector for all 
particles was in the same direction enabling a mid plane for 
the disc to be defined. Planet masses in the range 1-6 Mj 
were initiated in circular orbits with inclinations in the range 
io = [0; 80]° with respect to the initial disc mid plane. The 
semi-major axis of the planet was set to a = 5 AU. The 
particle distribution was chosen to model a surface density 
profile such that 



-1/2 



(3) 



Here Eo is a constant and R is the radial coordinate of a 
point in the midplane, here taken to occupy the radial do- 
main [0, Rout — S], with Rout = 5a and S = 0.4a. A taper 
was applied such that the surface density was set to decrease 
linearly to zero for R in the interval [Rout — S, Rout + S]. The 
total disc radial domain is hence [0, Rout + 5]. 
The disc mass is given by 



Md = 2n I nr)rdr = ^nEoRli] 



(4) 



which is used to determine Eq. The disc particles were set 
up in a state of pure Keplerian rotation velocity according 
to 



dr. 



(5) 



The radial velocities were set to be zero. Here $* is the grav- 
itational potential due to the central star. In the innermost 
region around the central star, the disc properties (e.g. Qs 
and accordingly Cs) are modified by the softening. 

In this paper, only discs for which self-gravity can be 
neglected are studied. For this to be possible with a locally 
isothermal equation of state, the Toomre parameter has to 
fulfil the requirement 



(6) 



For the case H/r^ — 0.05, the Toomre parameter can be 
expressed as 



Q 



0.05M. 

tIttE 



-^oMt^* ^-3/2 

15Md " 



(7) 



For a disc with any arbitrary outer radius Rout, the max- 
imum allowed disc mass M£),max, which fulfils the require- 
ment that Q{Rout) ~ 1.4 at its outer edge is 



Me 



M, 



l^Q{Rout) 21 



(8) 



Hence, for any arbitrarily chosen outer radius Rout, the total 
disc mass has to be set smaller than M«/21 ~ 0.05M«. For 
simulations shown in this paper, the disc mass is Md ~ 
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0.01 Mq unless stated otherwise which corresponds to 
3Md _ 0.03 M© 



So = 



4nRtil 20V57ra3/2 



(9) 



The disc mass surface density at the location of a planet 
with a = 5 AU is then 76 g cm~^. 

In addition to self-gravity, energy transfer and mass in- 
fall onto the disc are neglected. At the start of a simulation, 
the disc is allowed to evolve for 30 orbits before the planet is 
introduced. By doing this, we minimise the effects of initial 
transients in the disc arising from the setup. 

In order to examine the effect of changing the num- 
ber of SPH particles, we performed resolution studies with 
both coplanar planets and planets initiated on inclined or- 
bits. These are described in appendix[C] The number of SPH 
particles chosen for most of our simulations was 2x 10^. This 
enabled enough simulations to be run for ~ 500 orbits or 
500 X P(5 AU) = 500 x 11.18 yr = 5600 yr that a reasonable 
exploration of phase space could be made while retaining 
adequate accuracy. Finally we note that although for defi- 
niteness we have fixed length and time units to correspond 
to 5 AU, they can be scaled to other values with the usual 
corresponding scaling of the units of length and time. 



4 ORBITAL EVOLUTION OF A SINGLE 
MASSIVE PLANET 

4.1 General overview 

The evolution of a system composed of a central star, a 
planet and a gaseous disc is dependent on the relative masses 
of the different components. Low-mass planets produce rel- 
atively small amplitude density perturbations such that the 
global disc evolution is not significantly modified. In con- 
trast, a massive planet can influence the global evolution 
of a disc significantly in two ways. Firstly, material may 
be cleared from its neighbourhood by tidal interaction, cre- 
ating a deep and wide gap. Secondly, it is expected that 
massive planets on inclined orbits could modify the three- 
dimensional structure of the disc by producing large scale 
twists and warps. However, this has not yet been studied in 
detail. 

In the mass domain Alp £ [lMj;6Mj] considered here, 
a planet in a coplanar orbit is e xpected and been found 
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2006l l. In addition net 



to open a gap in t h e disc, (e.g. 
IVarniere etaP l2004l : ICrida et ail 
torques acting on the planet result in planet migration. 
When a planet opens a gap in its host disc, provided it is 
not too massive, it migrates inward on the viscous timescale 
of the disc (type II migration). 

When the orbit of the planet is significantly inclined 
with respect to the disc mid plane, the interaction with 
the disc is expected to differ from the coplanar case. Im- 
pulsive interactions occuring twice per orbit that can be de- 
scribed through an ap proach based on dynamical friction 
are expected. Recently, iReinI (|2012l ') made a local study of 
low-mass planets on highly inclined orbits and was able to 
confirm the dynamical friction viewpoint and show that the 
time-scale for realignment in such cases was long compared 
to the disc lifetime meaning that such inclined orbits should 
survive the presence of a disc. For massive planets, the in- 
teraction is expected to be stronger with the evolutionary 



time scales correspondingly reduced. In addition, the planet 
is more likely to affect the structure and evolution of the 
disc. 



4.2 Estimate of time scale for orbital evolution 

Here we estimate an approximate time scale for orbital evo- 
lution assuming that the interaction between the planet and 
the disc can be described as occurring through dynamical 
friction. It has been noted by a numb e r of authors (e.g. 
Papaloizou fc Larwoodl bOOCf: 'Reinl l2012l : [Teyssandier. et al.l 
2013) that planets on orbits with high eccentricity and/or 
high inclination pass through the disc with a high supersonic 
relative velocity such that the interaction can be approxi- 
mated as being with pressure-less ballistic fiuid particles. A 
dynamical friction calculation is then appropriate. The fric- 
tional force per unit mass is then (e.g. Rudermann fc Spieeell 
Il97ll : iRephaeli fc SalpetedlToSOl : IOstriked[l999l ) 



= i— 13 In ri/ra 



(10) 



where v^e; is the relative velocity between the planet and 
disc gas of density p. Here ri and r2 represent upper amd 
lower cut off length scales outside of which the analysis lead- 
ing to (|10[) that assumes a homogeneous medium becomes 
invalid. There is some uncertainty in estimating these that 
in turn leads to an uncertainty in the force of a factor of a 
few, and so we can only use (|10|l to make rough estimates. 
We consider a high inclination orbit and a thin disc such 
that the interaction may be viewed as impulsive and occur- 
ring twice per orbit. Each impulse produces a small velocity 
change given by 



Av = - 



4-!vGH'Ip\n{ri/r2) 



-dz, 



(11) 



where the disc mid plane is assumed to lie in a Cartesian 
{x,y) plane, the inclination of the orbit to the disc is i, the 
orbital speed is |v|, and the integration is taken through the 
vertical extent of the disc. For a thin disc this becomes 



Av : 



47rG^Mpln(r-i/r2)Evr 
|vre;P|v| sini 



(12) 



where E is the disc surface density. This in turn leads to a 
characteristic time for the orbit to change measured in units 
of the orbital period given by 



Td = |v,e;|/(2|Av|) = 



|v|* sin i sin'^(j/2) 
7rG'2Mpln(ri/r2)E' 



(13) 



where in order to obtain estimates, we have replaced \vrci\ 
by 2|v| sin(i/2). To proceed further we set r — a, jv|* — 
G'^M'i/a^ and E = 'iMoia/ Routf^'^ / {'^i^a'^)- Thus, we ob- 
tain 



A4",^sinisin^(i/2) _ AM'^ smism^ {i/2) f Rou. 



TrEMpflS ln(ri/r2) WIpMo ln(ri/r2) 



/ Rout \ 

\ a ) 



3/2 



(14) 



In order to estimate ri, we note that when the planet is 
in the disc, the vertical scale should be limited by H, but 
that r\ should be some what larger to account for more 
distant material near the disc midplane. Accordingly we take 
ri = 2H. For r2 we take the softening length which is around 
QAH at a = 5 AU. This is an estimate of the scale below 
which material should be regarded as being bound to the 
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Figure 1. The azimuthally averaged disc surface density pro- 
file (in units of Mj/(5AU)^) as a function of distance (in 
units of 5 AU) after 200 hundred orbits for planets with Mp = 
1,2,4, and 6 Mj. The curves plotted are such that the density 
minimum in the centre of the gap decreases with increasing planet 
mass. [All gap profiles were created using SPLASH tPriceii200l) ] 



planet. Taking Md/M* = 0.01, Mp/M. = 0.001, ri/r2 = 5, 
Rout /a = 5 and i = 7r/4, we obtain Td = 2 x 10* orbits, 
corresponding to 4.3 x 10^ yr for a = 5 AU. 



4.3 Lidov-Kozai effect 

In addition to dynamical friction, another mechanism that 
can affect the orbital evolution of planets on orbits that 
are highly in clined to the disc is the Lidov-Kozai effect 
l|Kozail 11963 '). This causes an exchange between inclina- 
tion and eccentricity for orbits above a critical inclina- 
tion. Orbits initially with high inclination can thus develop 
very high eccentricities. For the case of a stationary ax- 
isymmetric quadrupole pot ential, the critical i n clinat ion is 
estimated to be 3 9° w hile iTerguem fc Airnial (|2010l ) and 
iTevssandier et all (|2013l ') find values as small as ~ '20° for 
a planet interacting with an axisymmetric disc in two and 
three dimensional treatments respectively. In our case the 
disc is neither strictly axisymmetric nor steady. However, 
even though we find that it becomes warped and unsteady 
from the viewpoint of the planet, the effect remains rela- 
tively small so that to a first approximation, the planet may 
be considered as moving in an axisymmetric potential. How- 
ever, the component of angular momentum parallel to the 
symmetry axis will not be strictly conserved, as in that case, 
giving greater scope for generating inclination changes. 

It is important to note that in the axisymmetric case, in 
the secular the ory, a circular orb it at high inclination is on a 
separatrix (e.g. lFunk et al-lboilh . Thus it becomes an unsta- 
ble equilibrium point and chaotic motion is to be expected in 
its neighbourhood. This raises a potential difficulty in iden- 
tifying the cause of inclination changes associated with high 
inclination near circular orbits as being entirely due to fric- 
tional effects in simulation runs of limited length. However, 
we do not see significant eccentricity development, as would 
be expected if a Lidov-Kozai effect operated in such cases, 
which indicates that frictional effects are the main cause of 
the evolution. In addition the characteristic times for evo- 
lution can be compared to those expected from dynamical 
friction and they are found to be comparable. 
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Figure 2. Upper panel: gap formation for Mp = 1 Mj ini- 
tially with different equations of state, the azimuthally averaged 
surface density is shown after 200 orbits. The two upper curves 
corresponding to hp = 0.5 and hp = 0.6 are indistinguishable and 
the lowest curve is for the locally isothermal equation of state. 
Lower panel: mass growth starting from an initial planet mass 
Mp = 1 Mj for hp = 0.6 and locally isothermal equation of state. 
The upper curve corresponds to the locally isothermal equation 
of state. There is negligible accretion in the other case. 
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Figure 3. Upper panel: Semi-major axis as a function of time for 
coplanar planets of different masses (1, 2, 4, 6 Mj). Lower panel: 
cumulative mean migration rates as a function of time 
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Figure 4. Left: Relative inclination for planet masses (1, 2, 4, 6 Mj) as a function of time. These were initiated on prograde circular 
orbits with initial relative inclinations 10° , 20° , 30° , 40° , 60° and 80°. Right: Corresponding cumulative means of the relative inclination 
decay rates as a function of time. 



5 NUMERICAL RESULTS WITH A 
COPLANAR PLANET 

In order to understand the orbital evolution of massive plan- 
ets in the presence of a disc with its total angular momentum 
vector misaligned with that of the orbit, it is important to 
consider the extent of the occurrence of gap formation as 
this would expell material from the neighbourhood of the 
planets location and so is expected to affect the strength of 
the interaction of the planet with the disc. This is implied by 
the simple dynamical friction estimate for the force exerted 
between the planet and disc given by (|10|) . which indicates 
that this is proportional to the local disc density. Gap for- 



mation is most effective when the planet orbit and disc are 
coplanar as then the magnitude of the relative velocity be- 
tween the planet and disc material in its neighbourhood is 
the smallest. In addition, in contrast to the case of high or- 
bital inclination, the planet then spends most of its time 
close to disc material. Accordingly, we begin by studying 
the coplanar case. 

5.1 Gap formation in the coplanar case 

In this section, we consider the orbital evolution of planets 
with masses of 1, 2, 4 and 6 Mj initiated in coplanar cir- 
cular orbits with an initial semi-major axis of 5 AU. Fig. 
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[T] shows the azimuthaUy averaged disc surface density pro- 
file after 200 planet orbits with the modified equation of 
state. The development of a gap with a depth that increases 
with planet mass is clearly apparent. The viscous c riterion 
for gap formation given bv lLin fc Papaloizoul l|l993h is that 
Mp/M* > AQasshi, where for convenience, we have as- 
sumed that the gap extends beyond the region where the 
modified equation of state changes it from being locally 
is othermal and that dissipa tive processes are modelled by 
a lShakura fc SunvaevI (|l973l ) O-ss parameterization. For the 
simulations presented here, this has been estimated to be in 
the range 0.02 - 0.03 (see appendix E]) . Thus we should 
require Afp/M* > 0.002 - 0.003. This is consistent with 
Fig. [1] which indicates a gap depression of up to 66% for 
Mp/M. = 0.002 and a deep gap for Mp/M* > 0.004. The 
IMj case is associated with a partial gap associated with a 
33% depression of the azimuthaUy averaged surface density. 
In addition, we remark that lCrida et al.l (|200(i ') have studied 
the gap profile for different viscosities and planet masses. For 
Mp = 1 Mj, they find a gap depth of ~ 40% for hs = 0.05 
and ass = 0.025. We remark that the whole gap width we 
obtain, in this case being ab out 0.6 of the orbital radius, is 
very similar to that found bv lCrida et al. (2006 ). Our result s 
are also consistent with those of Ide Val-Borro et al.l (|2006l ) 
that indicated a somewhat less depleted gap found in SPH 
simulations as compared to grid based simulations. But note 
that the treatment of the equation of state and accretion of 
gas particles differs from ours. 

In order to illustrate the effects of the choice of the 
equation of state on the gap formation process and also the 
accretion of material by the planet, we illustrate a compar- 
ison between runs undertaken with the locally isothermal 
equation of state and the modified equation of state with 
hp — 0.5 and hp — 0.6 in Figure [21 

In the locally isothermal case, the planet mass increases 
to over 4 Mj on a time scale of about five hundred orbits. 
In this case, a high density of gas particles is able to ac- 
cumulate in the neighbourhood of the massive planet due 
to the small effective sound speed around the planet result- 
ing in relatively small pressure forces. Inside the Hill radius, 
SPH particles build up a hydrostatic equilibrium state with 
a very high mass density peak centred on the planet. The 
outer accretion radius is much smaller than the Hill radius 
in order to prevent any interference with the dynamics on 
that scale. Therefore, only particles that come very close to 
the planet can satisfy the conditions to be accreted. 

In contrast to the locally isothermal equation of state, 
the modified equation of state prevents gas particles from 
accumulating in the vicinity of the planet due to the sig- 
nificantly increased pressure forces occuring because of the 
increased sound speed. As a consequence, the Hill region 
around the planet does not contain a high density peak and 
accretion onto the planet is less significant. 



The results for h„ 



0.5 and h„ = 0.6 are almost iden- 



tical for the duration of the runs. Values hp < 0.5 were not 
used for 1 Mj since they result in regions around the planet 
where the modified Cs^mod is smaller than the Ca,iso as dis- 
cussed in appendixlX] It can be inferred from this study that 
the simulations shown in this paper that adopt the modified 
equation of state do not require implementation of the accre- 
tion algorithm. It is nonetheless retained for completeness 



and in order to ensure that any potential singularities at the 
locations of massive particles can be dealt with. 



5.2 Migration in the coplanar case 

In Fig. 131 the evolution of the semi-major axis a and the 
migration rate of a coplanar planet are studied. In order 
to eliminate any effects due to short timescale variations in 
defining migration rates, we evaluate a cumulative mean of 
da/dt as 



da 



da{t') _ a{0) - a(t) 
dt' ^ t 



(15) 



where a{t) is the semi-major axis at time t and a(0) is its 
value at t = 0. As expected, the migration is always found to 
be inwards. At the start of the simulations, no gap is present 
in the disc. During the gap formation phase, the migration 
rates are largest. At later times, the planets have opened 
gaps and the migration gradually slows down. The behaviour 
of the three largest masses is very similar. For Afp = 1 Mj, 
the gap is not as deep as in the other cases resulting in 
less of a decrease of the migration rate as compared to the 
other masses. At large times and for gap forming planets 
that are not too massive compared to the disc mass, inward 
migratio n is expected to occur o n the viscous timescale of 
the disc (|Lin fc Papaloizoul 1 19861 ). The stea dy state infiow 
velocity in a viscous disc is 1.5asshlrsns (e.g. |Pringlelll98ll ). 
This leads to an expected migration rate of Snasshl/ Porb, 
where Port is the local orbital period. For ass ~ 0.025, this 
corresponds to 2.5 x 10"* AU/yr at 5 AU. This is in good 
agreement with the rates shown at the latest time in Fig. O 
for the three larger masses. For Alp — 1 M,j, the migration 
rate is slightly slower at this time. This may be on account 
of the differing initial evolution, causing it to be at a larger 
radius and the presence of more material in the gap region. 



6 INTERACTIONS OF PLANETS WITH NON 
ZERO INCLINATION WITH THE DISC 

We now consider the evolution of planets initiated on cir- 
cular orbits with non zero relative inclination. The relative 
inclination, i, is defined through 



I = arccos 



Ip • L_D 



|lp| ■ 

Lu = ^ mi{ri X Vi). 



with 



(16) 
(17) 



Here Ip is the angular momentum vector of the planet and 
L_D is the total angular momentum vector of the disc, with 
the summation in p7|) being taken over all active gas parti- 
cles. In simulations starting with circular orbits, the plane- 
tary orbit develops only very small eccentricities over their 
run times. This is in contrast to the situation that occurs 
when orbits with high enough relative inclination are initi- 
ated with a modest non zero eccentricity (see Section (8}. 
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Figure 5. Gap formation for different planet masses and initial relative orbital inclinations: The azimuthally averaged surface density 
is shown after 200 orbits for, moving from left to right and top to bottom, Mp = 1,2,4 and 6 Mj, starting with different initial relative 
inclinations. 



6.1 Decay of the inclination for prograde initially 
circular orbits 

Fig. 3] shows the incUnation i of the planet orbit relative 
to the plane normal to the disc total angular momentum 
vector and the cumulative mean of its decay rate for planets 
of mass 1, 2, 4 and 6 Mj initiated on prograde circular orbits 
with initial relative inclinations 10° , 20° , 30° , 40° , 60° and 
80°. The cumulative mean of di/dt was adopted in order 
to iron out small time scale variations. This was defined 
through 

f^ifyj^dt'^m^, (18) 
dt t Jg dt' t ' ^ ^ 

with i{t) being the relative inclination at time t and j(0) its 
initial value. This approach loses meaning when the planet 
becomes almost coplanar with the disc. As can be seen from 
the left hand panels of Fig. |4l the latter readily happens for 
io < 20°. Cumulative means , defined as above would then 
approach i{0)/t, which is determined by the initial condi- 
tion. Accordingly decay rates are only shown for relative 
inclinations above 5° so avoiding this limit. 

In general, the cumulative mean of the inclination de- 
cay rate decreases in magnitude as the inclination increases. 
This qualitative feature is consistent with the simplistic 
estimate of the orbital evolution time scale based on dy- 
namical friction considerations given by HJ). This equa- 
tion suggests that this decay rate should be proportional 
to i/(sin i sin(i/2)'^). Thus comparing the early cumulative 
mean of the decay rates for i — 80° and i = 10° for 1 Mj 
in Fig. [l] we find a ratio of ~ 10"'^, with 1)14^ indicating a 
ratio of ~ 3 X 10~^. This ratio increases for larger masses 
on account of gap formation at low inclination. However, in 
all cases the ratio of the decay rates for i = 40° and i = 80° 
where gap formation is less relevant is ~ 16. In this case, 



(|14|l indicates a value of ~ 5. In view of the very rough na- 
ture of the estimate given by H14p the level of agreement is 
satisfactory. 

For io = 10° the results given in Fig. |4] indicate that 
Td ~ 10^ yr in the case of 1 Mj while for io = 80°, To ~ 
10^ yr initially. But note that the decay rate accelerates 
considerably as the inclination decreases so that for io = 40° , 
Td ~ 2 X 10^ yr . Again this is in reasonable agreement with 
the estimate given above derived from equation (|14p . 

For small to intermediate i and larger planet masses, the 
i-decay is strongly dependent on gap formation, i decreases 
at an accelerating rate towards smaller values as long as no 
gap is formed. As soon as i is sufficiently small to open a gap 
in the disc, the rate of inclination decrease slows down. At 
large inclinations, the inclination decay rate increases with 
planet mass so that for io = 80°, the decay rate is an order 
of magnitude faster for 6 Mj as compared to 1 Mj. 



6.2 The dependence of gap formation on orbital 
inclination for prograde orbits 

Because the relative velocity between the planet and disc 
material increases with orbital inclination, the interaction 
becomes weaker and gap formation is less likely. There can 
be a threshold inclination above which gap formation does 
not occur. This threshold inclination for gap opening in- 
creases with planet mass. Fig.[5]shows the azimuthally aver- 
aged disc surface mass density after 200 orbits for different 
planet masses and initial inclinations io. We remark that, 
as can be deduced from the results presented above, after 
200 orbits the relative inclination i will have decreased sig- 
nificantly from its initial value io in some cases. For very 
high inclinations, S is marginally perturbed for Mp = 6 Mj, 
the perturbation being weaker for the smaller masses as 
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Figure 6. Left: Evolution of the semi-major axis a for planets with Alp = 1 and 4 Mj for a range of initial relative orbital inclinations, 
Right: cumulative mean decay rates. 
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Figure 7. Orbital evolution of a planet with Mp = 4 Mj starting 
on retrograde circular orbits with different initial relative inclina- 
tions. In the upper panel, the semi-major axis is shown as function 
of time. Lower panel: relative inclination as function of time. 



expected. With decreasing io, the interaction between the 
planet and the disc becomes stronger until eventually the 
threshold below which a noticeable gap starts to form is 
passed. Thus for Mp = 1 Mj, only the io = 10° curve is able 
to open a noticeable partial gap. For Mp = 2 Mj,4 Mj and 
6 Mj, the threshold initial inclinations are io = 20°, 30° and 
40°, respectively. 



6.3 Evolution of the semi-major axis for finite 
relative orbital inclination 

Fig- El shows the evolution of the semi-major axes and cumu- 
lative means of the migration rates for planets with Mp = 1 
and 4 Mj initiated on circular orbits with different initial 
relative orbital inclinations. The cases Mp = 2 and 6 Mj 
show very similar behaviour, therefore will not be shown. 
The runs shown are identical to those illustrated in Fig. ID 
Thus corresponding relative inclinations can be read off. As 
in the coplanar case, the migration is always inwards and 
the orbits remain almost circular for the run times shown. 
For planets starting on inclined orbits, the migration rate is 
slower than for the coplanar case due to the reduced inter- 
action of the planet with the disc. When the relative orbital 
inclination approaches zero, it is found that the migration 
rate approaches the value appropriate to coplanar planets. 
For example, for Mp — 1 Mj and io — 20°, i decreases to 
< 5° after ~ 3000 years. During this time interval, the mi- 
gration rate approaches the value for a planet of the same 
mass with io = 0° . For Mp = 4 Mj and io — 30° , the relative 
inclination decreases to below 10° for t > 3000 yr. At this 
stage the migration rate approaches that for the coplanar 
case. 

6.4 High inclination retrograde orbits 

Finally in this section, we give examples of the orbital evolu- 
tion of a planet starting on a retrograde circular orbit with 
varying initial relative orbital inclination. In Fig. [T] we illus- 
trate the evolution of the semi-major axis and the relative 
orbital inclination as functions of time for a planet of mass 
Mp = 4 Mj. We remark that the direction of evolution of the 
inclination is always towards coplanarity {irei =0°), even 
when the orbit starts out with io ~ 180°. This is to be ex- 
pected from a frictional interaction between the planet and 
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Figure 8. Cross section of tlie disc as seen in the plane x 







z in units of 5 AU. The orientation of the coordinate system is 
such that the z axis corresponds to the direction of the angular 
momentum vector of the initial disc and the planet is initiated at 
y = with X > and z > 0. [All SPH cross section plots were 
created using SPLASH l |Pric6ll2007i) ] 




Figure 9. Cross section of the disc as seen in the piano x = 
after 200 orbits, for Mp = 1,2,4 and 6 Mj, and io = 40°. The 
coordinate system is as defined in Fig. [H] 



the disc which tends to communicate angular momentum in 
the direction of the disc's angular momentum vector to the 
orbit of the planet. In all of these cases, the orbit remains 
approximately circular. It is seen that the most rapid evolu- 
tion occurs for the largest inclinations for which the planet 
tends to become embedded in the disc and thus undergo a 
more sustained frictional interaction. 



7 THE RESPONSE OF THE DISC: 

INCLINATION CHANGES, WARPING AND 
PRECESSION 

As a response to a massive planet, a circumstellar disc 
can change its orientation significantly, developing a warped 
structure. So far, we have considered the inclination of the 
planetary orbits with respect to the plane for which the an- 




Figure 10. Cross section of the disc as seen in the plane x = 
after 200 orbits, for Mp = 1,2,4 and 6 Mj, and io = 80°. The 
coordinate system is as defined in Fig. \E\ 
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Figure 11. Evolution of total disc inclination for the first 5000 yr 
for different initial inclinations in the cases Mp = 1 and 6 Mj. 



gular momentum vector of the disc as a whole is parallel to 
the normal. However, warping can be present in the disc, 
dividing it into two parts with angular momentum vectors 
pointing in different directions and thus associated with dif- 
ferent inclinations. In Figure [S] we show a cross section con- 
sisting of the density distribution as seen in the plane x — 
for different planet masses, all of which started with an ini- 
tial inclination, io = 30°, after 200 orbits. This value of io 
was chosen as maximum warping is expected for an inclina- 
tion intermediate between 0° and 90°, there being none at 
these extremes. It is seen that a visible warped inner section 
of the disc is produced by the more massive planets while 
planets with Mp = 1 and 2 Mj do not succeed in creating a 
visibly warped disc. The same result is found for an initial 
inclination io — 40° as is illustrated in Fig. |9l The disc pro- 
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4 Mj and different initial relative inclinations for a 
~ 5000 years. The coordinate system is as defined in 
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Table 1. Characterising the disc warp; the first column gives 
the planet mass, the second column the initial relative orbital 
inclination, the third gives the relative orbital inclination after 200 
orbits and the fourth column gives the difference in the inclination 
angles of the inner and outer disc to the initial disc mid plane. 
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Figure 13. The magnitude of the retrograde precession angular 
velocity of the disc found from equation I I21I I as a function of 
time for planet masses with Mp = 1 and 6 Mj and various initial 
relative inclinations. 

in those cases. For the highest mass of 6 Mj, the total disc 
inclination only attains values up to 15°. 



file is shown for io — 80° for the different planet masses in 
Fig. IIUI In this case no warp is produced since the planets 
orbit is almost perpendicular to the disc mid plane. An in- 
dication of the degree of warping in the above cases can be 
obtained by measuring the inclination angles between the 
inner and outer discs as seen in these cross sections. It is 
noticeable that the inclination of the inner disc can be quite 
substantial in the direction of aligning with the inclination 
of the planetary orbit and that this effect is much less severe 
for the outer disc. Table [T] lists quantities characterising the 
disc warp after 200 orbits for three runs. As most of the 
inertia of the disc is in its outer parts, they dominate the in- 
clination of the disc defined by considering its total angular 
momentum. Fig. [TT]shows the evolution of the inclination of 
the disc, defined in this way, as a function of initial relative 
inclination of its orbit for the two extreme cases Mp = 1 
and Mp = 6 Mj. By definition the disc inclination is zero 
at the start of the simulations. The total disc inclination is 
determined from 



to 



D,z 



(19) 



In general, the disc inclination based on its total angular mo- 
mentum does not attain large values in any case because of 
the high disc inertia. For all the planet masses, the effect of 
the planet is strongest for the intermediate value io — 40° as 
expected. Simulations with io = 60° and io = 80° produce 
lower disc inclinations as would be expected for a rigid pla- 
nar disc because there would be reduced precessional torques 



7.1 Disc precession 

The gravitational interaction of a disc with a planet in an 
inclined orbit is expected to make the orbit precess about 
the total angular momentum vector of the system. Provided 
communication across the disc through either wave propa- 
gation or viscous diffusion occurs more rapidly than the pre- 
cession rate, the response of the disc would be expected to be 
to precess like a rigid body in a retrograde sense relative to 
the orbit (|Larwood et al.iri996l ) while simultaneously under- 
going a non disruptive warping. Due to angular momentum 
conservation of the composite system, this precession is also 
expected to be about the total angular momentum vector. 

In Figure [T^ we show the trajectory of the vector Lo.tot 
defined as 

hp X Ltot , . 

^OAot - -T^ — f 1 . (20) 

LiD X Litot\ 

Thus Li5,tot is a unit vector in the direction normal to the 
plane containing the disc angular momentum, Ld, and the 
total angular momentum vector, htot = \p + Lu, of the 
system. We focus on the case with Mp = 4 Mj, noting that 
other masses cases manifest similar behaviour. When there 
is rigid body precession, ho.tot rotates uniformly with the 
precession frequency about htot also in the retrograde sense. 
The magnitude of the angular velocity oiho.tot can be found 
from 



Ld tot X Vl 



(21) 



where v_l ~ dLD,tot/dt. This is shown as a function of time 
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Figure 14. Time-dependent evolution of the precession vector 
P(t) for Mp = 4 Mj and different initial inclinations io over 
~ 5000 yr. The coordinate system is as defined in Fig. [8] 



in Figure fTSl for the runs with planet masses Mp = 1, 6 Mj 
and initial relative orbital inclination. It is found to be al- 
most constant in all cases when a significant relative incli- 
nation is present that does not decrease significantly with 
time, which is a good indication of quasi rigid body preces- 
sion. It decreases in magnitude as the relative inclination 
increases which is expected as a reflection of the decrease 
in magnitude of the precessional torque acting between the 
planet and the disc as the relative incUnation increases. 

7.2 Precession of the line of nodes of inclined 
prograde planet orbits 

In addition to the disc precession, the planet orbit also pre- 
cesses. The orbital angular momentum vector should precess 
about the total angular momentum vector with the same 
angular velocity. Figure [14] shows the evolution of the pre- 
cession vector deflned as 



Ip X Jjtot 

Ip X Ltot 



(22) 



for the case Mp = 4 Mj and different initial inclinations for 
a time period of ~ 5000 yr. We remark that P is a unit 
vector perpendicular to the plane containing the orbital an- 
gular momentum vector and the total angular momentum 
vector. It is found to rotate with the same precession pre- 
cession frequency as the disc, as would be expected if the 
latter precessed like a rigid body . The same behaviour has 
also been found for all other planet mass cases studied. The 
precession is retrograde with the inclinations of the plane 
in which P evolves being dependent on the initial relative 
inclination of the planet. The higher the relative initial in- 
clination of the planet orbit, the higher the inclination of 
the plane of P. The angular velocity of orbital precession 
coincides exactly with that of the disc. 



8 EVOLUTION OF THE ECCENTRICITY AND 
THE INDICATION OF A LIDOV-KOZAI 
EFFECT AT HIGH INCLINATION 

So far, we have described simulations for which planets were 
initiated on circular orbits and in all cases the orbital eccen- 
tricity remained very small for the duration of the simula- 
tions. However, as the change in inclination of the major- 
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Figure 15. The orbital evolution for a Mp = 4 Mj planet starting 
on an eccentric orbit with e = 0.15 for initial relative orbital 
inclinations jq = 0° , 10° , 20° , 30° , 40° , 60° and 80° . 



ity of the disc to the initial mid plane and the correspond- 
ing warping distortions are relatively mild, the Lidov-Kozai 
mechanism applicable to orbits of high inclination with re- 
spect to the symmetry plane of an axisymmetric potential 
might be expected to operate (see the discussion in Section 
14.31 ) . As the part of a Lidov-Kozai cycle, where the orbit 
is nearly circular, corresponds to an extremum of the oscil- 
lation and being in the neighbourhood of a separatrix, and 
our runs were typically less than a precession cycle, devel- 
opment of an exchange between the inclination and eccen- 
tricity, characteristic of a Lidov-Kozai cycle, that could take 
place on longer time scales, may not have been observed. In 
order to clarify this issue, we initiated runs as before, but 
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starting with an initial orbital eccentricity of 0.15. In that 
case, exchange between the inclination and eccentricity is 
able to occur more rapidly such that it can be seen in our 
simulations for high masses and high inclinations. We show 
the results of such simulations for a planet with Mp — 4 M,j 
starting with different initial relative orbital inclinations in 
Fig. 1151 It is seen that the evolution of the semi-major axis 
proceeds in almost the same manner as for the cases that 
started with zero eccentricity. However, in this case we found 
that although for iq ^ 40°, the eccentricity ultimately de- 
cayed, for the cases of io ~ 60° and io = 80°, significant 
eccentricities > 0.4 developed and were increasing at the 
end of the simulations. This was accompanied by a decrease 
in relative inclination in those cases that was not observed 
in the cases initiated with circular orbits. 

Although these features are suggestive of the operation 
of a Lidov-Kozai like cycle for io >~ 40°, we remark that 
the observed decreases in inclination exceeded what would 
be expected from a conservation of the z component of an- 
gular momentum condition that would be applicable for an 
axisymmetric potential. This could be due to either the op- 
eration of dissipative decay and/or non axisymmetric effects. 
For io ^ 40° , the decay rates of the inclinations are about 
two to three times faster compared to those seen for simu- 
lations initiated with circular orbits. 

We also ran simulations for a 6 Mj planet starting with 
an initial eccentricity of 0.15 for different initial relative or- 
bital inclinations. Initial inclinations up to 60° show a simi- 
lar evolution to the 4 Mj case. For io = 80° , both the inclina- 
tion decrease and eccentricity increase occur more rapidly. 
After 5000 years, the relative inclination has decreased to 
40° while the eccentricity has reached 0.8. These features 
are also suggestive of the operation of a Lidov-Kozai like 
cycle. 

The characteristic precession period for initial relative 



orbital inclinations io ^ 60° can be extrapolated from Figs. 
1131 and [15] to be Pprec ~ 3.6 x 10* yr. This period might 
also be expected to be comparable to the period of Kozai 
oscillations. However, the warped structure of the precessing 
disc may play an important role (see Figs. [8] and [9| leading 
to modification of the standard description. 



8.1 Dependence on the disc mass 

So far, all simulations have been carried out for a standard 
disc mass of 0.01 Mq. It is expected that the interaction 
between planet and disc should increase with disc mass. 
Simple arguments based on dynamical friction that should 
be applicable to high inclination cases indicate that decay 
rates should scale with the disc mass. However, as the planet 
mass is in general not too different from the disc mass and 
a considerable amount of non-linearity and some distortion 
of the shape of the disc is present in our simulations, such a 
straightforward scaling may not be strictly valid. Figure [T6l 
shows the cumulative mean of relative inclination decay rate 
for Mp — 4 Mj initiated on a circular orbit with io = 60° for 
disc masses of 0.001,0.01, and 0.02 Mq. The first of these 
being 10 times smaller than the standard value and the last 
being twice as large. 

It is seen that the migration rates scale as the disc mass 
as expected from the simple dynamical friction argument. 
This is also the situation for the inclination decay rate when 
the larger of the two masses is considered. However, the 
smallest disc mass has a relative inclination decay rate that 
is 2 — 2.5 times faster than expected when compared to the 
standard mass. A closer scrutiny of these cases indicated 
that the shape of the disc differed in these two cases, with 
different orientations of the warp. This could result in rel- 
atively differing inclination changing torques that could be 
responsible for this result. 



9 CONCLUSIONS 

We have performed SPH-simulations in order to undertake 
a systematic study of the interaction of planets with masses 
in the range 1 — 6 Mj with initial inclinations in the range 
[0; 80]° with a circumstellar disc. For planets initiated on a 
coplanar circular orbit gap formation occurred. We found 
our gap profiles to be consistent with expectation s from the 
work of lLin fc Papaloizoul (| 19931 ) and|CridaeEaD (|2006l ) for 
grid based simulations for an adopted viscosity parameter 
ass = 0.025. Inward migration at a rate of ~ 2 x 10"* AU/yr 
for Mp =2, 4 and 6 Mj was found and this is in good agree- 
ment with previous studies and the theoretical expectation 
for type II migration. For Alp — 1 Mj, which was associ- 
ated with a partial gap, the migration rate was found to be 
slightly slower. 

We then considered the evolution of planets initiated 
on the full range of inclined circular orbits. For these sim- 
ulations, in contrast to those initiated with a finite orbital 
eccentricity, the orbits remained circular for the duration of 
the runs. It is possible that phenomena such as the Lidov- 
Kozai effect did not have time to operate in these cases. The 
inclination decay rate was found to decrease with increas- 
ing inclination and increasing planet mass. This is in full 
accordance with expectation from estimates based on the 
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operation of dynamical friction. For a disc mass of O.OIA/q, 
io = 80° and 6 Mj, the time to decay through 10° was 
found to be ~ 10^ yr while for 1 Mj, it is ~ 10*' yr. In 
the latter case, with the longest decay time, protoplanetary 
disc lifetimes are approached. For small to intermediate io, 
gap formation plays a crucial role in the determination of 
the inclination decay rates. Planets with larger masses are 
able to form gaps at higher relative inclinations. The re- 
duced amount of material near the planet then causes the 
frictional interaction to be reduced with a corresponding re- 
duction in the inclination decay rate. After 200 orbits, we 
find the threshold initial relative inclination below which 
gap formation starts to occur to be io = 10° for Mp = 1 Mj, 
io = 20° for Mp = 2 Mj, io = 30° for Mp = 4 Mj and 
io — 40° for Mp — 6 Mj, For planets on highly inclined or- 
bits, the migration rate was found to be reduced as compared 
to the coplanar case due to the reduced level of interaction 
with the disc. When the relative inclination approaches zero, 
the migration rate approaches that of the coplanar case. For 
the case of a 4 Mj planet initiated on a retrograde circular 
orbit and different initial inclinations, we found that the di- 
rection of evolution tends towards coplanarity (i = 0°) as 
would be expected from interaction with the disc through 
dynamical friction. The time scale for evolution decreased as 
io approached 180°. The most extreme case was found for 
io ~ 170° where the relative inclination decreased by ~ 25° 
in ~ 5000 years. 

We also found significant changes in the form of the disc 
that were produced by interaction with the larger planet 
masses. Visible warping occured in the inner parts. The dif- 
ference between the inclinations of the inner and outer part 
of the disc was found to attain up to 10 — 20° . Furthermore, 
we calculated the total disc orientation which can change 
by up to ~ 15° with respect to its original mid plane. As 
expected, we also found retrograde precession of both the 
the total disc angular momentum vector and the planet's 
angular momentum vector about the conserved angular mo- 
mentum vector of the system. This occurred at an almost 
constant angular velocity over simulation times, with a mag- 
nitude that decreased with increasing relative inclination. 

In addition to simulations starting with zero eccen- 
tricity, we initiated a 4 Mj planet on inclined orbits with 
e = 0.15. While for io < 40°, the eccentricity decayed, the 
cases with io = 60 and 80° showed an eccentricity increasing 
for the duration of the simulations. The eccentricity increase 
was accompanied by a relatively rapid decrease in relative 
inclination. These characteristics indicate the possibility of 
a Lidov-Kozai like effect operating in these cases. 

In this paper, we have made a first attempt to quantify 
the possible influences of massive planets with a range of dif- 
ferent initial orbital inclinations on a circumstellar disc. Our 
results are reasonably consistent with theoretical estimates 
of evolution time scales from dynamical friction, although 
the simulated behaviour of the disc turns out to be complex 
in some cases and was not taken into account. In contrast 
to the situation with low-mass planets, whose influence on 
a disc is negligible, the results shown in this paper suggest 
that the interaction of a massive planet with a disc can lead 
to a strong distortion. Furthermore the stronger frictional 
interaction is likely to lead to coplanarity within disc life- 
times in most cases. Thus highly inclined orbits are only 



likely to survive if they are formed after the disc has mostly 
dispersed. 
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APPENDIX A: CHOICE OF ASPECT RATIO IN 
THE CIRCUMPLANETARY DISC FOR THE 
MODIFIED EQUATION OF STATE 

The main motivation of IPeplinski et al] (|2008l ) is to intro- 
duce a modification to the locally isothermal equation of 
state in order to prevent a large amount of gas from being 
accreted by massive planets. In the close vicinity of massive 
planets, the sound speed is enhanced, as would be expected 
for accreted optically thick material in hydrostatic equilib- 
rium. This feature is not described by a locally isothermal 
equation of state. Thus the modified sound speed Cs,mod is 
designed to be the sum of the unmodified soundspeed Cs,iso 
and an additional positive quantity whose value increases 
towards the planet and tends to zero far from it. 

Given the motivation indicated above, Cs.mod should be 
everywhere larger than Cs ,iso- [Peplinski et al. (2008:) have 
shown through various tests with planets with mass ex- 
ceeding IMj, that hp ^ 0.4 is appropriate in order to pre- 
vent large amounts of gas from being accreted. Our studies 
have shown that for a given planet mass, Mp, there exists 
a lower limit hpjim which must be exceeded in order to 
guarantee that cs.mod ^ Cs.iso globally. For hp ^ /ip,iim, 
Cs,mod ^ Cs,iso somcwhcrc and the modification does not 
comply with physical requirements. In Table lAll some values 
of hp,iirn are listed. The lower bound hp^um decreases with 
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Figure Al. Modified sound speed and locally isothermal sound 
speed for Mp = 0.1 Mj, hp = 0.8 (top); Mp = 1 Mj, hp = 0.4 
(middle); Mp = 1 Mj, hp = 0.6 (bottom). 
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Table Al. Lower bounds hp n^ for different planet masses 

increasing planet mass. For 1 Mj, we suggest hp > 0.49 for 
the above reasons. Fig. lAll illustrates two examples where 
Ca,mod < Cs,iso and & third that achieves a successful modi- 
fication of the equation of state. 

APPENDIX B: RING SPREADING TEST 

The spreading of a viscous gas ring is often used to calibrate 
the effective viscosity (e.g. ISpeith fc Riffert|[l999l ). The an- 
alytic solution describing an axisymmetric thin disc of gas 
moving around a central point mass Mc under gravity and 
viscous stresses is considered. The vertical scale height, H, 
of the disc is assumed to be negligible compared to the ra- 
dial length scale and the viscous time scale governing the 
radial inflow of meiss is assumed to be much larger than the 
dynamical time scale. In the thin-disc approximation, it is 
possible to integrate the governing equations over the height 
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Figure Bl. Comparison between the analytical and numerical 
solutions with different values for the artificial viscosity parameter 
a (see Section 12.21 1 as indicated. 



of the disc and consider a two dimensional flow. Adopting 
polar coordinates {R, ip), there are two velocity components 
{vr,v^), and the mass distribution is described by the sur- 
face density 



pdz 



(Bl) 



where p is the mass density. Assuming that the radial in- 
flow of mass is caused entirely by viscous stresses, the sur- 
face density c an be shown to satisfy the diffusion equation 
l|Pringlelll98"ll ) 



9E _ 3__9_ 

dt ^ RdR 
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(B2) 



If the initial surface density distribution is chosen to be pro- 
portional to a (5 function corresponding to a ring of mass 
M localised at an initial radius iioi then the initial surface 
density distribution T,o{R) is given by 

M 

MR)-^J{R-Ro) (B3) 

The analytic solution of (|B2|) gives E at a later time as 
M 1 , /2a;\ / I 



T.{R,t) = 



with vr = 



Rot 



C-/-3/4(^)//l/4(^ 



(B4) 
(B5) 



Here 7 are the modified Bessel functions, x — R/Ro, 
T = V2vQt/R^ is the time in units of the viscous ti me scale 
fln/(' 12i^o) and vq is the kinematic viscosity (see iPringld 
Il98lb . We adopt the values Rq = 1 Rq, M = 10"^° M©, 
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Figure CI. Resolution study of gap formation (upper panel) and 
corresponding migration (lower panel) for a coplanar planet with 
Mp = 4 Mj and a disc of Mo = 0.01 Mq. Curves corresponding 
to runs with different numbers of particles as indicated. 



Mc = 1 M0, I/O = 1.5 • lO^"* cmVs (see also ISpeith fc Kiev! 
(|2003h ). Since E(t = 0) is a 5-function, the test is started by 
initiating the state variables to match conditions at r = 0.02 
and then the system is evolved forward s in time. The corre- 
spon ding dimensionless ass parameter (jShakura fc SunvaevI 
Il973t ) corresponding to this optimised value of vq is 



ass = 



CsH 



1.5 ■ 10" cmVs 

0.052 ^/gm|1; 



= 0.0196 ^ 0.02 



(B6) 



Results from our test calculations with 2 x 10^ SPH parti- 
clcs, arc shown in Figure IBll In contrast to ISpeith fc Kiev! 
(2003), we didn't find any indications for spiral arm forma- 
tion. In this context we note that according to these authors, 
the spiral features are associated with instabilities driven by 
shear viscosity and that they used a version of SPH with a 
treatment of viscosity that differs from the one we used. 



APPENDIX C: RESOLUTION STUDY 

In order to investigate the sensitivity of the simulations to 
the number of SPH particles used, we have performed a 
number of studies with differing numbers of particles for a 
planet with Mp = 4 Mj initiated on both circular orbits, and 
orbits with an initial eccentricity, around a central star with 
M« = IM0 and interacting with a disc with Md ~ O.OIMq. 
We considered the coplanar case, io — 0, and the case where 
io = 60°. 

Runs were performed with, 10^,2 x 10^, and 4 x 10^ 
SPH particles. We remark that as shown by ring spreading 
tests, a disc modelled with SPH particles has a non zero 
effective shear viscosity that is expected to decrease weakly 
as the number of particles is increased. Thus one wants to 
verify that a fluid disc with an effective shear viscosity is 
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Figure C2. Resolution study of relative inclination evolution 
(upper panel) and inward migration (lower panel) for a planet 
with Mp = 4 Mj and iq = 60°. The disc mass was Mo = 
0.01 Mq. Curves corresponding to runs with different numbers 
of particles as indicated. 



appropriately modelled for the number of particles chosen. 
For all these test studies, the disc was allowed to evolve for 
30 orbits without the planet in order to allow any initial 
inhomogenities and fluctuations to be smoothed out. The 
planet was then suddenly introduced after 30 orbits. The 
sudden introduction produces transients that rapidly decay 
and so produce no lasting noticeable effects. 

In Fig. ICll we show the results for initiation on a circu- 
lar orbit with io = 0. Gap profiles and the evolution of the 
semi-major axis are plotted for the three runs with differing 
numbers of SPH particles. For both gap profiles and the evo- 
lution of the semi-major axis, the different resolution levels 
show good agreement over the run time of the simulations. 
There is a hint of a decreasing shear viscosity with increas- 
ing numbers of particles, as the semi-major axis decreases 
slightly less, but this effect is small. 

Fig. IC2I shows the evolution of the relative inclination 
and the semi-major axis for the cases when the planet is 
initiated on a circular orbit with io = 60°. The results of 
simulations carried out with different numbers of particles 
are also seen to be in good agreement. In this case the evolu- 
tion is expected to be driven by dynamical friction and thus 
expected to be insensitive to disc viscosity. The evolution of 
the semi-major axis is indistinguishable in these cases. The 
correspondence for the evolution of the relative inclination 
is not quite as precise. However, note that the evolution of 
the relative inclination is particularly sensitive because the 
initial orbit is likely to be close to a separatrix and its initial 
evolution thus expected to be sensitive to small perturba- 
tions (see discussion in Section 

We have also run a resolution test for a planet starting 
with io = 60° and an initial eccentricity e = 0.15. This study 
was performed to check the results described in Section [S] 
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Figure C3. Top two panels: as for Fig. IC2l but starting with an 
initial eccentricity of 0.15. Bottom panel: The evolution of the 
eccentricity is shown. 



The evolution is again expected to be controlled by dynam- 
ical friction and insensitive to disc viscosity. The inclination 
and semi-major axis evolution are plotted in Fig. IC3I The 
results of simulations carried out with different numbers of 
particles are seen to be in good agreement. The runs under- 
taken with different numbers of particles thus indicate good 
convergence and that a choice of 2 x 10^ SPH particles for 
simulations of the type presented here is a reasonable one. 

This paper has been typeset from a T^gX/ ffl^jX file prepared 
by the author. 



